/* ======Trend Inflation and the Nature of Structural Breaks in the New Keynesian Phillips Curve ======

This program estimates the bivariate unobserved components model for the NKPC in inflation and real output
with two structural breaks. The variability of shocks to the trend component of inflation is constrained
to be constant.
 
Written by: Pym Manopimoke

Model:

Pi_t = Pi_t^* + k_{S_t} Sum_{j=0}^{inf} beta^j E_{t-1}(x_{t+j}) + z_t
Pi_t^* = Pi_{t-1}^* + v_t, v_t ~iid.N(0,sig^2_v)
z_t = psi_{S_t} z_{t-1} + epsilon_t, epsilon_t|S_t ~iid.N(0,sig^2_{epsilon,S_t})
y_t = tau_t + x_t
tau_t = mu*D_{1t} + tau_{t-1} + eta_t, eta_t|D_{1t} ~iid.N(0,sig^2{eta,D_{2t}})
x_t = phi_1 x_{t-1} + phi_2 x_{t-2} + u_t, u_t|D_{2t}~iid.N(0,sig^2_{u,D_{2t}))

u_t and epsilon_t are correlated with correlation coefficient rho.

S_t = {0,1,2}
Prob[ S(t)=0 | S(t-1)=0] = p00
Prob[ S(t)=1 | S(t-1)=1] = p11
Prob[ S(t)=2 | S(t-1)=2] = 1

mu_{D_{1t}} = (1-D_{1t})*mu_0 + D_{1t}*mu_1
D_{1t} = {0 if t<= 1973Q1, 1 otherwise}
sig^2_{u,D_{2t}} = sig^2_{u,0}(1-D_{2t}) + sig^2_{u,1}D_{2t}
D_{2t} = {0 if t<= 1984Q3, 1 otherwise}

========================================================================================================*/

new;
library optmum,PGRAPH;
pqgwin auto;
cls;
 
load data[250,11]=data.prn;  
data=data[2:250,.];
t=rows(data);
@Set the sample to 1952:1-2007:3@
startsmpl=21;
endsmpl =243; 
  
qtr= data[startsmpl:endsmpl,1];  

/*Column 2 for GDP chained-weighted inflation, 7 for CPI inflation, 8 for PCE inflation, 
9 for GDP Deflator*/  
inf=data[startsmpl:endsmpl,8];  @PCE inflation@
yy=ln(data[startsmpl:endsmpl,3])*100; 
xx=data[startsmpl:endsmpl,5];  @CBO output gap@
t=rows(inf~yy~xx);

@Define dummy variable series D_1t and D_2t to account for productivity slowdown and the Great Moderation@
PS=85; @1973.1@
GM = 130; @1984:2@

dmat = ones(GM,1)|zeros(t-GM,1); 
pmat=ones(PS,1)|zeros(t-PS,1);

@================= Initialize Global Variables============@

tvp_no=5; 
start= 11;
t1=75;t2=125;
   
prmtr_in=  4.3287975|
           3.6169770|
           0.01|
           0.01|
           1.0257639|
           0.46744571|
           0.93990821|
           0.74457509|
           0.33874759|
           0.74178097|
           1.5487628|
           0.84717542|
           -0.076618725|
           1.5183078|
          -5.4242774|
          -0.069540940|
           8.7701265|
           0.0049308014|
           0.021651444|
           0.0043142585|
           0.013086457;
    
{xout,fout,gout,cout}=optmum(&lik_fcn,PRMTR_IN);
PRM_FNL=TRANS(xout);             @ Estimated  coefficients, constrained@
xout;
hout0=hessp(&lik_fcn,xout);
hout0;
hout=inv(hout0);
grdn_fnl=gradfd(&TRANS,xout);
Hsn_fnl=grdn_fnl*hout*grdn_fnl';
SD_fnl =sqrt(diag(Hsn_fnl)); @Standard errors of the estimated coefficients@
 
{B_TT,PR_TT,PR_TL,stdprederr,xfor,confint} =FILTER(XOUT);
smooth3=SMOOTH(pr_tt,pr_tl,prm_fnl[1:2,1]); @Pr[S_t|Y_T], Smoothed probabilities@

"== FINAL OUTPUT ===============================";
"Likelihood value is ";; -fout;
"Estimated parameters (s.e.) tstat are:"; 
prm_fnl~sd_fnl~prm_fnl./sd_fnl;
"===============================================";

@Produces Figure A3 in the unpublished appendix@
_pdate=0;
_pltype = {5 6 3};
_plegctl={4 4 6 5};
_plegstr="Regime1\000Regime2\000Regime3";
title("Smoothed probabilities of regimes");
title("");
XY(qtr[start:t],SMOOTH3); 

@Produces Figure A4 in the unpublished appendix@
_pltype = {6 3 5 5};
_plegctl={2 4 5.5 5.3};
_plegstr="Trend Inflation\000Actual Inflation\000Lower and Upper\00090\% Confidence Bands";
title("Actual inflation and trend inflation");
title("");
xy(qtr[start:t],b_tt[.,4]~inf[start:t]~confint);

@Produces Figure A5 in the unpublished appendix@
_pltype = {6 3};
_plegctl={2 4 6 5.7};
_plegstr="Inflation Cycle\000Inflation Gap";
title("Inflation cycle and inflation gap");
title("");
XY(qtr[start:t],xfor~inf[start:t]-b_tt[.,4]);

@Produces Figure A6 in the unpublished appendix@
_pltype = {6 3};
_plegctl={2 4 4.8 5.7};
_plegstr="Bivariate UC Model Output Gap\000CBO Output Gap";
title("Output gaps");
title("");
xy(qtr[start:t],b_tt[.,2]~xx[start:t]);
  
output file=UC1.out reset;
 
output off;
end;

@ end OF MAIN PROGRAM @
@========================================================================@
 

PROC LIK_FCN(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
sigetaGM0,sigetaGM1,mu0,mu1,pt,mu,sigeta,
xfor,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2, 
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22; 


EXTERNAL PROC TRANS, V_PROB;
 
           PRMTR=TRANS(PRMTR1);

           P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
           P11=PRMTR[2,1];    @Pr[St=1/St-1=1]@
           P22=1;    @Pr[St=2/St-1=2]@
           sigetaGM0=prmtr[3,1];
           sigetaGM1=prmtr[4,1];
           siguGM0=prmtr[5,1];
           siguGM1=prmtr[6,1];
           mu0=prmtr[7,1];
           mu1=prmtr[8,1];
           sigv=prmtr[9,1];  
           sigep=prmtr[10:12,1];  
           rho=prmtr[13,1];
           phi1=prmtr[14,1];
           phi2=prmtr[15,1];
           psi=prmtr[16:18,1];    
           k=prmtr[19:21,1];

           A=1~0;
           H = 0~0~0~1~1|1~1~0~0~0;
           
           B_LL0 = yy[1]|0|0|inf[1]|0;      @initial guess for coefficients@
           P_LL0 = eye(5)*10000;      @uncertainty associated with initial guess@
           B_LL1=yy[t1]|0|0|inf[t1]|0; 
           P_LL1 = eye(5)*10000;       
           B_LL2=yy[t2]|0|0|inf[t2]|0; 
           P_LL2 = eye(5)*10000;      
                
           subF = phi1~phi2| 1~0; 
           subQ= siguGM0^2~0|0~0;
           vecp = inv(eye(4)-subF.*.subf)*vec(subq);
           vv = reshape(vecp',2,2);
           p_ll0[2:3,2:3] = vv;
           p_ll1[2:3,2:3] = vv;
           p_ll2[2:3,2:3] = vv;
           p_ll0[5,5] = sigep[1,1]^2/(1-psi[1]^2);
           p_ll1[5,5] = sigep[2,1]^2/(1-psi[2]^2);
           p_ll2[5,5] = sigep[3,1]^2/(1-psi[3]^2);

           PROB__0 = 1; @Unconditional probability for St=0@
           PROB__1 = 0; @Unconditional probability for St=1@
           PROB__2 = 0; @Unconditional probability for St=2@
 

/*  START ITERATION */
/*===================================================================== */

           LIKV=0.0;
           J_ITER = 1;
           do until J_ITER>T;

                yt=inf[j_iter]|yy[j_iter];
                pt=pmat[j_iter]; 
                mu=mu0*pt+mu1*(1-pt);
                dmt = dmat[j_iter];
                sigu=siguGM0*dmt + siguGM1*(1-dmt);
                sigeta=sigetaGM0*dmt + sigetaGM1*(1-dmt);

                Q0=sigeta^2~0~0~0~0|
                   0~sigu^2~0~0~rho*sigu*sigep[1,1]|
                   0~0~0~0~0|
                   0~0~0~sigv^2~0|
                   0~rho*sigu*sigep[1,1]~0~0~sigep[1,1]^2;
                     
                Q1=sigeta^2~0~0~0~0|
                   0~sigu^2~0~0~rho*sigu*sigep[2,1]|
                   0~0~0~0~0|
                   0~0~0~sigv^2~0|
                   0~rho*sigu*sigep[2,1]~0~0~sigep[2,1]^2;
                     
                Q2=sigeta^2~0~0~0~0|
                   0~sigu^2~0~0~rho*sigu*sigep[3,1]|
                   0~0~0~0~0|
                   0~0~0~sigv^2~0|
                   0~rho*sigu*sigep[3,1]~0~0~sigep[3,1]^2;

                F0 = 1~0~0~0~0|
                     0~phi1~phi2~0~0|
                     0~1~0~0~0|
                     0~0~0~1~0|
                     0~0~0~0~psi[1];
                F1 = 1~0~0~0~0|
                      0~phi1~phi2~0~0|
                      0~1~0~0~0|
                      0~0~0~1~0|
                      0~0~0~0~psi[2];
                F2 = 1~0~0~0~0|
                      0~phi1~phi2~0~0|
                      0~1~0~0~0|
                      0~0~0~1~0|
                      0~0~0~0~psi[3];

                FF= phi1~phi2|1~0;
                xfor = k[1]*A*(inv(eye(2)-FF));           
                ome0 = 0~xfor~0~0| 0~0~0~0~0;
                xfor = k[2]*A*(inv(eye(2)-FF)); 
                ome1 = 0~xfor~0~0| 0~0~0~0~0;
                xfor = k[3]*A*(inv(eye(2)-FF)); 
                ome2 = 0~xfor~0~0| 0~0~0~0~0;

                /* KALMAN FILTER */
               /*==============================================================*/
                 @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
                 B_TL00 = F0 * B_LL0+(mu|0|0|0|0);  @WHEN S_{t-1}=0, S_{t}=0@
                 B_TL01 = F1 * B_LL0+(mu|0|0|0|0);  @WHEN S_{t-1}=0, S_{t}=1@
                 B_TL02 = F2 * B_LL0+(mu|0|0|0|0);  @WHEN S_{t-1}=0, S_{t}=2@
                 B_TL10 = F0 * B_LL1+(mu|0|0|0|0);  @WHEN S_{t-1}=1, S_{t}=0@
                 B_TL11 = F1 * B_LL1+(mu|0|0|0|0);  @WHEN S_{t-1}=1, S_{t}=1@
                 B_TL12 = F2 * B_LL1+(mu|0|0|0|0);  @WHEN S_{t-1}=1, S_{t}=2@
                 B_TL20 = F0 * B_LL2+(mu|0|0|0|0);  @WHEN S_{t-1}=2, S_{t}=0@
                 B_TL21 = F1 * B_LL2+(mu|0|0|0|0);  @WHEN S_{t-1}=2, S_{t}=1@
                 B_TL22 = F2 * B_LL2+(mu|0|0|0|0);  @WHEN S_{t-1}=2, S_{t}=2@
    
                 @ P_tl(i,j),  S(t-1)=i, S(t)=j@
                 P_TL00 = F0* P_LL0 * F0' + Q0;
                 P_TL01 = F1* P_LL0 * F1' + Q1;
                 P_TL02 = F2* P_LL0 * F2' + Q2;
                 P_TL10 = F0* P_LL1 * F0' + Q0;
                 P_TL11 = F1* P_LL1 * F1' + Q1;
                 P_TL12 = F2* P_LL1 * F2' + Q2;
                 P_TL20 = F0* P_LL2 * F0' + Q0;
                 P_TL21 = F1* P_LL2 * F1' + Q1;
                 P_TL22 = F2* P_LL2 * F2' + Q2;
 
                 @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
                 F_CAST00= yt- H * B_TL00 - ome0*B_TL00;
                 F_CAST01= yt- H * B_TL01 - ome1*B_TL01;
                 F_CAST02= yt- H * B_TL02 - ome2*B_TL02;
                 F_CAST10= yt- H * B_TL10 - ome0*B_TL10;
                 F_CAST11= yt- H * B_TL11 - ome1*B_TL11 ;
                 F_CAST12= yt- H * B_TL12 - ome2*B_TL12 ;
                 F_CAST20= yt- H * B_TL20 - ome0*B_TL20;
                 F_CAST21= yt- H * B_TL21 - ome1*B_TL21;
                 F_CAST22= yt- H * B_TL22 - ome2*B_TL22 ;
 
                 @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
                 SS00= H* P_TL00 * H' ; 
                 SS01= H* P_TL01 * H' ;
                 SS02= H* P_TL02 * H' ;
                 SS10= H* P_TL10 * H' ;
                 SS11= H* P_TL11 * H' ;
                 SS12= H* P_TL12 * H' ;
                 SS20= H* P_TL20 * H' ;
                 SS21= H* P_TL21 * H' ;
                 SS22= H* P_TL22 * H' ;
                 
                 /*Updating*/
                 B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
                 B_TT01 = B_TL01 + (P_TL01 * H' *inv( SS01)) * F_CAST01;
                 B_TT02 = B_TL02 + (P_TL02 * H' *inv( SS02)) * F_CAST02;
                 B_TT10 = B_TL10 + (P_TL10 * H' *inv( SS10)) * F_CAST10;
                 B_TT11 = B_TL11 + (P_TL11 * H' *inv( SS11)) * F_CAST11;
                 B_TT12 = B_TL12 + (P_TL12 * H' *inv( SS12)) * F_CAST12;
                 B_TT20 = B_TL20 + (P_TL20 * H' *inv(SS20)) * F_CAST20;
                 B_TT21 = B_TL21 + (P_TL21 * H' *inv(SS21)) * F_CAST21;
                 B_TT22 = B_TL22 + (P_TL22 * H' *inv(SS22)) * F_CAST22;

                 P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
                 P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
                 P_TT02 = (EYE(tvp_no) - (P_TL02 * H' *inv(SS02) * H )) * P_TL02;
                 P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
                 P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
                 P_TT12 = (EYE(tvp_no) - (P_TL12 * H' *inv(SS12) * H )) * P_TL12;
                 P_TT20 = (EYE(tvp_no) - (P_TL20 * H' *inv(SS20) * H )) * P_TL20;
                 P_TT21 = (EYE(tvp_no) - (P_TL21 * H' *inv(SS21) * H )) * P_TL21;
                 P_TT22 = (EYE(tvp_no) - (P_TL22 * H' *inv(SS22) * H )) * P_TL22;
       
                 /* HAMILTON FILTER */      

                 PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
                 PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
                 PR_VL02=0;
                 PR_VL10=0;
                 PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
                 PR_VL12=V_PROB(F_CAST12,SS12)*(1-p11)*PROB__1;
                 PR_VL20=0;
                 PR_VL21=0;
                 PR_VL22=V_PROB(F_CAST22,SS22)* p22 *PROB__2;
                
                 @CONDITIONAL DENSITY TIMES WEIGHT@

                 PR_VAL=PR_VL00+PR_VL01+PR_VL02+PR_VL10+PR_VL11+PR_VL12+PR_VL20+PR_VL21+PR_VL22;
                 
                 @Pr[St,St-1/Yt]@
                 PRO_00=PR_VL00/PR_VAL;       
                 PRO_01=PR_VL01/PR_VAL;
                 PRO_02=PR_VL02/PR_VAL;
                 PRO_10=PR_VL10/PR_VAL;
                 PRO_11=PR_VL11/PR_VAL;
                 PRO_12=PR_VL12/PR_VAL;
                 PRO_20=PR_VL20/PR_VAL;
                 PRO_21=PR_VL21/PR_VAL;
                 PRO_22=PR_VL22/PR_VAL;
                 PROB__0=PRO_00+PRO_10+PRO_20;       @Pr[St=0/Yt]@
                 PROB__1=PRO_01+PRO_11+PRO_21;       @Pr[St=1/Yt]@
                 PROB__2=PRO_02+PRO_12+PRO_22;       @Pr[St=2/Yt]@


       /* COLLAPSING TERMS */
       /*==============================================================*/

              if(prob__0!=0);
              B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 + PRO_20*B_TT20)/PROB__0;
              endif;
              if(prob__1!=0);
              B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 + PRO_21*B_TT21)/PROB__1;
              endif;
              if(prob__2!=0);
              B_LL2=(PRO_02*B_TT02 + PRO_12*B_TT12 + PRO_22*B_TT22)/PROB__2;
              endif;
              if(prob__0!=0);
              P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                     PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)')+
                     PRO_20*(P_TT20+(B_LL0-B_TT20)*(B_LL0-B_TT20)'))/PROB__0;
              endif;
              if(prob__1!=0);
              P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                     PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)')+
                     PRO_21*(P_TT21+(B_LL1-B_TT21)*(B_LL1-B_TT21)'))/PROB__1;
              endif;
              if(prob__2!=0);
              P_LL2=(PRO_02*(P_TT02+(B_LL2-B_TT02)*(B_LL2-B_TT02)')+
                     PRO_12*(P_TT12+(B_LL2-B_TT12)*(B_LL2-B_TT12)')+
                     PRO_22*(P_TT22+(B_LL2-B_TT22)*(B_LL2-B_TT22)'))/PROB__2;
              endif;
              if J_ITER < START; goto skip; endif;
 
              LIKV = LIKV -LN(PR_VAL);
              skip:
       J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
  
RETP(LIKV);
ENDP;

 
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
proc TRANS(prmtr1); @ constraining values of reg. coeff.@
local prm_new,d0,d1;

prm_new=prmtr1; 

prm_new[1:2,1]=EXP(prmtr1[1:2,1]) ./ (1 + EXP(prmtr1[1:2,1])); 
prm_new[16:18] =prmtr1[16:18,1]./(1+abs(prmtr1[16:18,1]));
prm_new[13] =prmtr1[13,1]./(1+abs(prmtr1[13,1]));
prm_new[3:6,1] = sqrt(prmtr1[3:6,1]^2);        
prm_new[9:12,1] = sqrt(prmtr1[9:12,1]^2);
d0 =prmtr1[14,1]./(1+abs(prmtr1[14,1]));
d1 =(1-abs(d0))*prmtr1[15,1]./(1+abs(prmtr1[15,1]))+ abs(d0)-d0^2;
prm_new[14,1]=2*d0;
prm_new[15,1]=-1*(d0^2 + d1) ;
  
retp(prm_new);
endp;

@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@
PROC V_PROB(EV, HE);      @ CALCULATES    Pr[Yt/St,Yt-1]@
local VAL;
         
VAL=(1/SQRT(2*PI*det(HE)))*EXP(-0.5*EV'inv(HE)*EV);

RETP(VAL);
ENDP;
@>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>@


proc (6) = FILTER(PRMTR1);

local PRMTR,p00,p11, p22,siguGM0,siguGM1,sigu,sigv,sigep,rho,psi,PHI1,PHI2,k,
sigetaGM0,sigetaGM1,mu0,mu1,pt,mu,sigeta,
xfor0,xfor1,xfor2,dmt,ome,ome0,ome1,ome2,subf,subq,H,vv,vecp,FF,A,yt,confint,lowerCI,upperCI,stderr,stdprederr,xforecast,
PROB__0,PROB__1,PROB__2,LIKV,J_ITER,Q0,Q1,Q2,F0,F1,F2,B_M,PR_TTM,PR_TLM,
pr_tl00,pr_tl01, pr_tl02, pr_tl10,pr_tl11,pr_tl12,pr_tl20,pr_tl21,pr_tl22,
B_TL00,B_TL01,B_TL02,B_TL10,B_TL11,B_TL12,B_TL21,B_TL20,B_TL22,     
P_TL00,P_TL01,P_TL02,P_TL10,P_TL11,P_TL12,P_TL20,P_TL21, P_TL22, 
F_CAST00,F_CAST01,F_CAST02,F_CAST10,F_CAST11,F_CAST12,F_CAST20,F_CAST21,F_CAST22,
SS00,SS01,SS02,SS10,SS11,SS12,SS20,SS21,SS22,
B_LL0,B_LL1,B_LL2,P_LL0,P_LL1,P_LL2,pr_ll0_,pr_ll1_pr_ll2_, 
B_TT00,B_TT01,B_TT02,B_TT10,B_TT11,B_TT12,B_TT20,B_TT21,B_TT22,
P_TT00,P_TT01,P_TT02,P_TT10,P_TT11,P_TT12,P_TT20,P_TT21,P_TT22,
PR_VL00,PR_VL01,PR_VL02,PR_VL10,PR_VL11,PR_VL12,PR_VL20,PR_VL21,PR_VL22,
PR_VAL,PRO_00,PRO_01,PRO_02,PRO_10,PRO_11,PRO_12,PRO_20,PRO_21,PRO_22;  

EXTERNAL PROC TRANS, V_PROB;

           PRMTR=TRANS(PRMTR1);
 
           P00=PRMTR[1,1];    @Pr[St=0/St-1=0]@
           P11=PRMTR[2,1];    @Pr[St=1/St-1=1]@
           P22=1;    @Pr[St=2/St-1=2]@
           sigetaGM0=prmtr[3,1];
           sigetaGM1=prmtr[4,1];
           siguGM0=prmtr[5,1];
           siguGM1=prmtr[6,1];
           mu0=prmtr[7,1];
           mu1=prmtr[8];
           sigv=prmtr[9,1];  
           sigep=prmtr[10:12,1];  
           rho=prmtr[13,1];
           phi1=prmtr[14,1];
           phi2=prmtr[15,1];
           psi=prmtr[16:18,1];    
           k=prmtr[19:21,1];

           A=1~0;
           H = 0~0~0~1~1|1~1~0~0~0;
           
           B_LL0 = yy[1]|0|0|inf[1]|0;      @initial guess for coefficients@
           P_LL0 = eye(5)*10000;      @uncertainty associated with initial guess@
           B_LL1=yy[t1]|0|0|inf[t1]|0; 
           P_LL1 = eye(5)*10000;       
           B_LL2=yy[t2]|0|0|inf[t2]|0; 
           P_LL2 = eye(5)*10000;      
                
           subF = phi1~phi2| 1~0; 
           subQ= siguGM0^2~0|0~0;
           vecp = inv(eye(4)-subF.*.subf)*vec(subq);
           vv = reshape(vecp',2,2);
           p_ll0[2:3,2:3] = vv;
           p_ll1[2:3,2:3] = vv;
           p_ll2[2:3,2:3] = vv;
           p_ll0[5,5] = sigep[1,1]^2/(1-psi[1]^2);
           p_ll1[5,5] = sigep[2,1]^2/(1-psi[2]^2);
           p_ll2[5,5] = sigep[3,1]^2/(1-psi[3]^2);

           PROB__0 = 1; @Unconditional probability for St=0@
           PROB__1 = 0; @Unconditional probability for St=1@
           PROB__2 = 0; @Unconditional probability for St=2@
 
           B_M=ZEROS(T,tvp_no);  @will store (1,1)element of B_TT@
           PR_TTM=ZEROS(T,3);  @will store Pr[S_t=0|Y_t]@
           PR_TLM=ZEROS(T,3);  @will store Pr[S_t=0|Y_{t-1}]@
           
           stdprederr=zeros(t,1); @standardized prediction error for diagnostic tests@
           confint=zeros(t,2);  @90% confidence interval bands for trend inflation estimates@ 
           xforecast=zeros(t,1); @stores the infinite sum term in the cyclical component@
          

/*  START ITERATION */
/*===================================================================== */

       LIKV=0.0;
       J_ITER = 1;
       do until J_ITER>T;
                 yt=inf[j_iter]|yy[j_iter];    
                 pt=pmat[j_iter]; 
                 mu=mu0*pt+mu1*(1-pt);
                 dmt = dmat[j_iter];
                 sigu=siguGM0*dmt + siguGM1*(1-dmt);
                 sigeta=sigetaGM0*dmt + sigetaGM1*(1-dmt);

                 pr_tl00 = PROB__0*p00;         
                 pr_tl01 = PROB__0*(1-p00);
                 pr_tl02 = 0;
                 pr_tl10 = 0;
                 pr_tl11 = PROB__1*p11;
                 pr_tl12 = PROB__1*(1-p11);
                 pr_tl20 = 0;
                 pr_tl21 = 0;
                 pr_tl22 = p22 *PROB__2;
                 PR_TLM[J_ITER,.]= (pr_tl00 +pr_tl10+pr_tl20)
                         ~(pr_tl01 +pr_tl11+pr_tl21)
                         ~(pr_tl02 +pr_tl12+pr_tl22); 

                 Q0=   sigeta^2~0~0~0~0|
                       0~sigu^2~0~0~rho*sigu*sigep[1,1]|
                       0~0~0~0~0|
                       0~0~0~sigv^2~0|
                       0~rho*sigu*sigep[1,1]~0~0~sigep[1,1]^2;
                         
                 Q1=   sigeta^2~0~0~0~0|
                       0~sigu^2~0~0~rho*sigu*sigep[2,1]|
                       0~0~0~0~0|
                       0~0~0~sigv^2~0|
                       0~rho*sigu*sigep[2,1]~0~0~sigep[2,1]^2;
                         
                 Q2=   sigeta^2~0~0~0~0|
                       0~sigu^2~0~0~rho*sigu*sigep[3,1]|
                       0~0~0~0~0|
                       0~0~0~sigv^2~0|
                       0~rho*sigu*sigep[3,1]~0~0~sigep[3,1]^2;
    
                 F0 =  1~0~0~0~0|
                       0~phi1~phi2~0~0|
                       0~1~0~0~0|
                       0~0~0~1~0|
                       0~0~0~0~psi[1];
                 F1 =  1~0~0~0~0|
                       0~phi1~phi2~0~0|
                       0~1~0~0~0|
                       0~0~0~1~0|
                       0~0~0~0~psi[2];
                 F2 =  1~0~0~0~0|
                       0~phi1~phi2~0~0|
                       0~1~0~0~0|
                       0~0~0~1~0|
                       0~0~0~0~psi[3];
    
                 FF= phi1~phi2|1~0;
                 xfor0 = k[1]*A*(inv(eye(2)-FF));           
                 ome0 = 0~xfor0~0~0| 0~0~0~0~0;
                 xfor1 = k[2]*A*(inv(eye(2)-FF)); 
                 ome1 = 0~xfor1~0~0| 0~0~0~0~0;
                 xfor2 = k[3] *A*(inv(eye(2)-FF)); 
                 ome2 = 0~xfor2~0~0| 0~0~0~0~0;


       /* KALMAN FILTER */
       /*==============================================================*/

                 @ Beta_tl(i,j),  S(t-1)=i, S(t)=j @
                 B_TL00 = F0 * B_LL0+(mu|0|0|0|0);  @WHEN S_{t-1}=0, S_{t}=0@
                 B_TL01 = F1 * B_LL0+(mu|0|0|0|0);  @WHEN S_{t-1}=0, S_{t}=1@
                 B_TL02 = F2 * B_LL0+(mu|0|0|0|0);  @WHEN S_{t-1}=0, S_{t}=2@
                 B_TL10 = F0 * B_LL1+(mu|0|0|0|0);  @WHEN S_{t-1}=1, S_{t}=0@
                 B_TL11 = F1 * B_LL1+(mu|0|0|0|0);  @WHEN S_{t-1}=1, S_{t}=1@
                 B_TL12 = F2 * B_LL1+(mu|0|0|0|0);  @WHEN S_{t-1}=1, S_{t}=2@
                 B_TL20 = F0 * B_LL2+(mu|0|0|0|0);  @WHEN S_{t-1}=2, S_{t}=0@
                 B_TL21 = F1 * B_LL2+(mu|0|0|0|0);  @WHEN S_{t-1}=2, S_{t}=1@
                 B_TL22 = F2 * B_LL2+(mu|0|0|0|0);  @WHEN S_{t-1}=2, S_{t}=2@
    
                 @ P_tl(i,j),  S(t-1)=i, S(t)=j@
                 P_TL00 = F0* P_LL0 * F0' + Q0;
                 P_TL01 = F1* P_LL0 * F1' + Q1;
                 P_TL02 = F2* P_LL0 * F2' + Q2;
                 P_TL10 = F0* P_LL1 * F0' + Q0;
                 P_TL11 = F1* P_LL1 * F1' + Q1;
                 P_TL12 = F2* P_LL1 * F2' + Q2;
                 P_TL20 = F0* P_LL2 * F0' + Q0;
                 P_TL21 = F1* P_LL2 * F1' + Q1;
                 P_TL22 = F2* P_LL2 * F2' + Q2;
 
                 @ Forecast error: F_cast_tl(i,j),  S(t-1)=i, S(t)=j@
                 F_CAST00= yt- H * B_TL00 - ome0*B_TL00;
                 F_CAST01= yt- H * B_TL01 - ome1*B_TL01;
                 F_CAST02= yt- H * B_TL02 - ome2*B_TL02;
                 F_CAST10= yt- H * B_TL10 - ome0*B_TL10;
                 F_CAST11= yt- H * B_TL11 - ome1*B_TL11 ;
                 F_CAST12= yt- H * B_TL12 - ome2*B_TL12 ;
                 F_CAST20= yt- H * B_TL20 - ome0*B_TL20;
                 F_CAST21= yt- H * B_TL21 - ome1*B_TL21;
                 F_CAST22= yt- H * B_TL22 - ome2*B_TL22 ;
 
                 @VARIANCE OF FORECAST ERROR: SS_tl(i,j),  S(t-1)=i, S(t)=j@
                 SS00= H* P_TL00 * H' ; 
                 SS01= H* P_TL01 * H' ;
                 SS02= H* P_TL02 * H' ;
                 SS10= H* P_TL10 * H' ;
                 SS11= H* P_TL11 * H' ;
                 SS12= H* P_TL12 * H' ;
                 SS20= H* P_TL20 * H' ;
                 SS21= H* P_TL21 * H' ;
                 SS22= H* P_TL22 * H' ;
                 
                 /*Updating*/
                 B_TT00 = B_TL00 + (P_TL00 * H' *inv(SS00)) * F_CAST00;
                 B_TT01 = B_TL01 + (P_TL01 * H' *inv( SS01)) * F_CAST01;
                 B_TT02 = B_TL02 + (P_TL02 * H' *inv( SS02)) * F_CAST02;
                 B_TT10 = B_TL10 + (P_TL10 * H' *inv( SS10)) * F_CAST10;
                 B_TT11 = B_TL11 + (P_TL11 * H' *inv( SS11)) * F_CAST11;
                 B_TT12 = B_TL12 + (P_TL12 * H' *inv( SS12)) * F_CAST12;
                 B_TT20 = B_TL20 + (P_TL20 * H' *inv(SS20)) * F_CAST20;
                 B_TT21 = B_TL21 + (P_TL21 * H' *inv(SS21)) * F_CAST21;
                 B_TT22 = B_TL22 + (P_TL22 * H' *inv(SS22)) * F_CAST22;

                 P_TT00 = (EYE(tvp_no) - (P_TL00 * H' *inv(SS00) * H )) * P_TL00;
                 P_TT01 = (EYE(tvp_no) - (P_TL01 * H' *inv(SS01) * H )) * P_TL01;
                 P_TT02 = (EYE(tvp_no) - (P_TL02 * H' *inv(SS02) * H )) * P_TL02;
                 P_TT10 = (EYE(tvp_no) - (P_TL10 * H' *inv(SS10) * H ) )* P_TL10;
                 P_TT11 = (EYE(tvp_no) - (P_TL11 * H' *inv(SS11) * H )) * P_TL11;
                 P_TT12 = (EYE(tvp_no) - (P_TL12 * H' *inv(SS12) * H )) * P_TL12;
                 P_TT20 = (EYE(tvp_no) - (P_TL20 * H' *inv(SS20) * H )) * P_TL20;
                 P_TT21 = (EYE(tvp_no) - (P_TL21 * H' *inv(SS21) * H )) * P_TL21;
                 P_TT22 = (EYE(tvp_no) - (P_TL22 * H' *inv(SS22) * H )) * P_TL22;

                /* HAMILTON FILTER */      

                 PR_VL00=V_PROB(F_CAST00,SS00)*p00 *PROB__0;    @Pr[St,Yt/Yt-1]@
                 PR_VL01=V_PROB(F_CAST01,SS01)*(1-p00)*PROB__0;
                 PR_VL02=0;
                 PR_VL10=0;
                 PR_VL11=V_PROB(F_CAST11,SS11)*p11 *PROB__1;
                 PR_VL12=V_PROB(F_CAST12,SS12)*(1-p11)*PROB__1;
                 PR_VL20=0;
                 PR_VL21=0;
                 PR_VL22=V_PROB(F_CAST22,SS22)* p22 *PROB__2;
                
                 @CONDITIONAL DENSITY TIMES WEIGHT@

                 PR_VAL=PR_VL00+PR_VL01+PR_VL02+PR_VL10+PR_VL11+PR_VL12+PR_VL20+PR_VL21+PR_VL22;
                 
                 @Pr[St,St-1/Yt]@
                 PRO_00=PR_VL00/PR_VAL;       
                 PRO_01=PR_VL01/PR_VAL;
                 PRO_02=PR_VL02/PR_VAL;
                 PRO_10=PR_VL10/PR_VAL;
                 PRO_11=PR_VL11/PR_VAL;
                 PRO_12=PR_VL12/PR_VAL;
                 PRO_20=PR_VL20/PR_VAL;
                 PRO_21=PR_VL21/PR_VAL;
                 PRO_22=PR_VL22/PR_VAL;

                 PROB__0=PRO_00+PRO_10+PRO_20;       @Pr[St=0/Yt]@
                 PROB__1=PRO_01+PRO_11+PRO_21;       @Pr[St=1/Yt]@
                 PROB__2=PRO_02+PRO_12+PRO_22;       @Pr[St=2/Yt]@
                 PR_TTM[J_ITER,.]=PROB__0~prob__1~prob__2;

                 /* COLLAPSING TERMS */
      
                 if(prob__0!=0);
                 B_LL0=(PRO_00*B_TT00 + PRO_10*B_TT10 + PRO_20*B_TT20)/PROB__0;
                 endif;
                 if(prob__1!=0);
                 B_LL1=(PRO_01*B_TT01 + PRO_11*B_TT11 + PRO_21*B_TT21)/PROB__1;
                 endif;
                 if(prob__2!=0);
                 B_LL2=(PRO_02*B_TT02 + PRO_12*B_TT12 + PRO_22*B_TT22)/PROB__2;
                 endif;
                 if(prob__0!=0);
                 P_LL0=(PRO_00*(P_TT00+(B_LL0-B_TT00)*(B_LL0-B_TT00)')+
                 PRO_10*(P_TT10+(B_LL0-B_TT10)*(B_LL0-B_TT10)')+
                 PRO_20*(P_TT20+(B_LL0-B_TT20)*(B_LL0-B_TT20)'))/PROB__0;
                 endif;
                 if(prob__1!=0);
                 P_LL1=(PRO_01*(P_TT01+(B_LL1-B_TT01)*(B_LL1-B_TT01)')+
                 PRO_11*(P_TT11+(B_LL1-B_TT11)*(B_LL1-B_TT11)')+
                 PRO_21*(P_TT21+(B_LL1-B_TT21)*(B_LL1-B_TT21)'))/PROB__1;
                 endif;
                 if(prob__2!=0);
                 P_LL2=(PRO_02*(P_TT02+(B_LL2-B_TT02)*(B_LL2-B_TT02)')+
                 PRO_12*(P_TT12+(B_LL2-B_TT12)*(B_LL2-B_TT12)')+
                 PRO_22*(P_TT22+(B_LL2-B_TT22)*(B_LL2-B_TT22)'))/PROB__2;
                 endif;
                 
                 @Estimates of the unobserved components@
                 B_TT=    B_TT00*PRO_00+B_TT01*PRO_01+B_TT02*PRO_02+
                          B_TT10*PRO_10+B_TT11*PRO_11+B_TT12*PRO_12+
                          B_TT20*PRO_20+B_TT21*PRO_21+B_TT22*PRO_22;
                 B_M[J_ITER,.]=B_TT';


                @90% Confidence bands for trend inflation estimates@
                lowerCI=(B_TT00[4,1]-1.645*sqrt(P_TT00[4,4]))*PR_tl00+
                       (B_TT01[4,1]-1.645*sqrt(P_TT01[4,4]))*PR_tl01+
                       (B_TT02[4,1]-1.645*sqrt(P_TT02[4,4]))*PR_tl02+
                       (B_TT10[4,1]-1.645*sqrt(P_TT10[4,4]))*PR_tl10+
                       (B_TT11[4,1]-1.645*sqrt(P_TT11[4,4]))*PR_tl11+
                       (B_TT12[4,1]-1.645*sqrt(P_TT12[4,4]))*PR_tl12+
                       (B_TT20[4,1]-1.645*sqrt(P_TT20[4,4]))*PR_tl20+
                       (B_TT21[4,1]-1.645*sqrt(P_TT21[4,4]))*PR_tl21+
                       (B_TT22[4,1]-1.645*sqrt(P_TT22[4,4]))*PR_tl22;
                upperCI=(B_TT00[4,1]+1.645*sqrt(P_TT00[4,4]))*PR_tl00+
                       (B_TT01[4,1]+1.645*sqrt(P_TT01[4,4]))*PR_tl01+
                       (B_TT02[4,1]+1.645*sqrt(P_TT02[4,4]))*PR_tl02+
                       (B_TT10[4,1]+1.645*sqrt(P_TT10[4,4]))*PR_tl10+
                       (B_TT11[4,1]+1.645*sqrt(P_TT11[4,4]))*PR_tl11+
                       (B_TT12[4,1]+1.645*sqrt(P_TT12[4,4]))*PR_tl12+
                       (B_TT20[4,1]+1.645*sqrt(P_TT20[4,4]))*PR_tl20+
                       (B_TT21[4,1]+1.645*sqrt(P_TT21[4,4]))*PR_tl21+
                       (B_TT22[4,1]+1.645*sqrt(P_TT22[4,4]))*PR_tl22;
                confint[j_iter,.] = lowerCI~upperCI;

                @Standardized prediction error@
                stderr= (f_cast00[1,1]/sqrt(ss00[1,1]))*pr_tl00+ (f_cast01[1,1]/sqrt(ss01[1,1]))*pr_tl01+
                       (f_cast02[1,1]/sqrt(ss02[1,1]))*pr_tl02+ (f_cast10[1,1]/sqrt(ss10[1,1]))*pr_tl10+
                       (f_cast11[1,1]/sqrt(ss11[1,1]))*pr_tl11+ (f_cast12[1,1]/sqrt(ss12[1,1]))*pr_tl12+ 
                       (f_cast20[1,1]/sqrt(ss20[1,1]))*pr_tl20+ (f_cast21[1,1]/sqrt(ss21[1,1]))*pr_tl21+
                       (f_cast22[1,1]/sqrt(ss22[1,1]))*pr_tl22;
                STDPREDERR[j_iter,.] = stderr;
               
                @infinite sum term in the cyclical component for inflation@
                xforecast[j_iter,.]=(xfor0 *B_TL00[2:3,1])*pr_tl00 +
                                    (xfor1 *B_TL01[2:3,1])*pr_tl01 +
                                    (xfor2 *B_TL02[2:3,1])*pr_tl02 +
                                    (xfor0 *B_TL10[2:3,1])*pr_tl10 +
                                    (xfor1 *B_TL11[2:3,1])*pr_tl11 +
                                    (xfor2 *B_TL12[2:3,1])*pr_tl12 + 
                                    (xfor0 *B_TL20[2:3,1])*pr_tl20 +
                                    (xfor1 *B_TL21[2:3,1])*pr_tl21 +
                                    (xfor2 *B_TL22[2:3,1])*pr_tl22;

              IF J_ITER < START; goto skip; endif;
              LIKV = LIKV -LN(PR_VAL);
              skip:
              J_ITER = J_ITER+1;
       ENDO;

/*  END ITERATION    */
/*===================*/

retp(B_M[START:T,.],PR_TTM[START:T,.],PR_TLM[START:T,.],stdprederr[start:t,.],xforecast[start:t,.],confint[start:t,.]);
ENDP;


@==========================================================@

@Kim's smoothing@
proc SMOOTH(prtt,prtl,prmtr0);

local   pr_sm1,pr_sm2,pr_sm3, j_iter, TT,prmtr, pr_tr, p11,p12,p13,p21,p22,p23,p31,p32,p33,
        pr_tt1,pr_tt2,pr_tt3,pr_tl1,pr_tl2,pr_tl3,pr_sm11,pr_sm12,pr_sm13,pr_sm21,pr_sm22,pr_sm23,pr_sm31,pr_sm32,pr_sm33;


     TT=ROWS(PRTT);
     prmtr=prmtr0;
 
     p11=prmtr[1,1]; 
     p12=1-p11; 
     p13=0; 
     p21=0; 
     p22=prmtr[2,1]; 
     p23=1-p22;
     p31=0; 
     p32=0; 
     p33=1;
 
     pr_tt1=prtt[.,1];
     pr_tt2=prtt[.,2];
     pr_tt3=prtt[.,3];
     pr_tl1=prtl[.,1];
     pr_tl2=prtl[.,2];
     pr_tl3=prtl[.,3];
           
     pr_sm1=pr_tt1;      @ pr_sm1 will contain Pr[S_t=1|Y_T]@
     pr_sm2=pr_tt2;      @ pr_sm2 will contain Pr[S_t=2|Y_T]@
     pr_sm3=pr_tt3;      @ pr_sm3 will contain Pr[S_t=3|Y_T]@
          

     j_iter=TT-1;
     do until j_iter < 1;
        if pr_tl1[j_iter+1,1]<0.000001; pr_tl1[j_iter+1,1]=0.0000001; endif;
        if pr_tl2[j_iter+1,1]<0.000001; pr_tl2[j_iter+1,1]=0.0000001; endif;
        if pr_tl3[j_iter+1,1]<0.000001; pr_tl3[j_iter+1,1]=0.0000001; endif;
         
        pr_sm11=pr_sm1[j_iter+1,1]*p11*  pr_tt1[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm12=pr_sm2[j_iter+1,1]*p12*  pr_tt1[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm13=pr_sm3[j_iter+1,1]*p13*  pr_tt1[j_iter,1]/ pr_tl3[j_iter+1,1];
 
        pr_sm21=pr_sm1[j_iter+1,1]*p21*  pr_tt2[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm22=pr_sm2[j_iter+1,1]*p22*  pr_tt2[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm23=pr_sm3[j_iter+1,1]*p23*  pr_tt2[j_iter,1]/ pr_tl3[j_iter+1,1];
      
        pr_sm31=pr_sm1[j_iter+1,1]*p31*  pr_tt3[j_iter,1]/ pr_tl1[j_iter+1,1];
        pr_sm32=pr_sm2[j_iter+1,1]*p32*  pr_tt3[j_iter,1]/ pr_tl2[j_iter+1,1];
        pr_sm33=pr_sm3[j_iter+1,1]*p33*  pr_tt3[j_iter,1]/ pr_tl3[j_iter+1,1];
       
        pr_sm1[j_iter,1]=pr_sm11+pr_sm12+pr_sm13;
        pr_sm2[j_iter,1]=pr_sm21+pr_sm22+pr_sm23;
        pr_sm3[j_iter,1]=pr_sm31+pr_sm32+pr_sm33;
    
     j_iter=j_iter -1;
     endo;
 
retp(pr_sm1~pr_sm2~pr_sm3);  
endp;

/*=================================================================================================================*/
 